Loading the libraries and data

I’ll first load all relevant libraries.

library (osmdata)
library (opentripplanner)
library (tidyverse)
library (sf)
library (ggthemes)
library (ggspatial)

Load locations

I’ll load in the locations for Boston public schools, and then filter:

  1. For only public schools in Dorchester
  2. For only public elementary (K-3) schools in Dorchester
BPS_Geojson <- st_read ("http://bostonopendata-boston.opendata.arcgis.com/datasets/1d9509a8b2fd485d9ad471ba2fdb1f90_0.geojson?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D")

Dorchester_BPS <- BPS_Geojson %>%
  filter (CITY == "Dorchester")

Dorchester_EL_BPS <- BPS_Geojson %>%
  filter (CITY == "Dorchester", SCH_TYPE == "ES")

Get street data

I’ll get street data for Dorchester, Boston, MA.

MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

dot_street_query <- opq(bbox = 'Dorchester MA USA')%>%
  add_osm_feature (key = 'highway')

dot_street_query %>%
  osmdata_xml (file = 'OTP/Graphs/Default/dot_streets.osm')

dot_street_features <- dot_street_query %>%
  osmdata_sf()

dot_streets <- dot_street_features$osm_lines %>%
  st_transform(crs=MA_state_plane)

Plot Dorchester streets

ggplot(dot_streets) + 
  geom_sf () +
  theme_map ()

Set up OpenTripPlanner

path_otp <- otp_dl_jar("OTP")
path_data <-file.path (getwd(), "OTP")
path_otp <- paste (path_data, "otp.jar", sep ="/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024) 
## 2020-10-05 14:14:05 Basic checks completed, building graph, this may take a few minutes
## The graph will be saved to C:/Users/melja/OneDrive/Desktop/GSD/VIS/melmiller-vis/OTP
## 2020-10-05 14:14:35 Graph built

Launch OpenTripPlanner

otp_setup(otp = path_otp, dir = path_data, memory =1024)

Connect to server

otpcon <- otp_connect()
## Router http://localhost:8080/otp/routers/default exists

Create isochrones

I’ll first create isochrones for areas within a five-minute bike and five-minute walk from all Boston Public Schools in Dorchester.

iso_5min_walk <-
    otp_isochrone(otpcon = otpcon, fromPlace = Dorchester_BPS, 
                mode = "WALK", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")

iso_5min_bike <- 
  otp_isochrone(otpcon = otpcon, fromPlace = Dorchester_BPS, 
                mode = "BICYCLE", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "bike")

iso_all_modes <- rbind(iso_5min_bike, iso_5min_walk)
right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = Dorchester_BPS) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("by bike", "by foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors") +
  theme(legend.position= "right")
## Loading required namespace: raster

From this point forward I’ll focus on just elementary schools (K-3) in Dorchester. Note that I’ve chosen to not include K-8 and K-12 schools in Dorchester.

iso_5min_walk <-
    otp_isochrone(otpcon = otpcon, fromPlace = Dorchester_EL_BPS, 
                mode = "WALK", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")

iso_5min_bike <- 
  otp_isochrone(otpcon = otpcon, fromPlace = Dorchester_EL_BPS, 
                mode = "BICYCLE", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "bike")

iso_all_modes_EL <- rbind(iso_5min_bike, iso_5min_walk)

otp_stop()
## [1] "SUCCESS: The process \"java.exe\" with PID 12260 has been terminated."
right_side <- st_bbox(iso_all_modes_EL)$xmax
left_side  <- st_bbox(iso_all_modes_EL)$xmin
top_side <- st_bbox(iso_all_modes_EL)$ymax
bottom_side <- st_bbox(iso_all_modes_EL)$ymin

ggplot(iso_all_modes_EL) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = Dorchester_EL_BPS) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("by bike", "by foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors") +
  theme(legend.position= "right")

Thunderforestlandscape basemap

While I think this map is too busy, I think it’s interesting to see the isochrones contextualized with the major transit stops and parks/open space.

right_side <- st_bbox(iso_all_modes_EL)$xmax
left_side  <- st_bbox(iso_all_modes_EL)$xmin
top_side <- st_bbox(iso_all_modes_EL)$ymax
bottom_side <- st_bbox(iso_all_modes_EL)$ymin

ggplot(iso_all_modes_EL) +
  annotation_map_tile(zoomin = 0, type="thunderforestlandscape", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = Dorchester_EL_BPS) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_discrete(name = "Area that is reachable within 5 minutes",
                       labels = c("by bike", "by foot"),
                       type = c("white", "black")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors") +
  theme(legend.position= "right")

Cartolight basemap

This map provides an understanding of major roads and neighborhood squares (e.g., Fields Corner). I think this map does the best job so far at helping visualize where K-3 elementary schools are, along with their accessibility (by bike and foot).

right_side <- st_bbox(iso_all_modes_EL)$xmax
left_side  <- st_bbox(iso_all_modes_EL)$xmin
top_side <- st_bbox(iso_all_modes_EL)$ymax
bottom_side <- st_bbox(iso_all_modes_EL)$ymin

ggplot(iso_all_modes_EL) +
  annotation_map_tile(zoomin = 0, type="cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = Dorchester_EL_BPS) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("by bike", "by foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors") +
  theme(legend.position= "right")

Hillshade basemap

This last map uses shading to show hills, or the lack of, in Dorchester.

right_side <- st_bbox(iso_all_modes_EL)$xmax
left_side  <- st_bbox(iso_all_modes_EL)$xmin
top_side <- st_bbox(iso_all_modes_EL)$ymax
bottom_side <- st_bbox(iso_all_modes_EL)$ymin

ggplot(iso_all_modes_EL) +
  annotation_map_tile(zoomin = 0, type="hillshade", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = Dorchester_EL_BPS) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("by bike", "by foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors") +
  theme(legend.position= "right")

Calculate and compare isochrone areas

I’ll calculate the area of the biking and walking isochrones to visualize the relationship between the two.

iso_areas <- iso_all_modes_EL %>%
  mutate(area = st_area(iso_all_modes_EL)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(bike))) +
  geom_point() +
  scale_x_continuous(name = 
            "Area within a five-minute walking distance\nof a public elementary (K-3) school\n(square km)",
            breaks = breaks <- seq(10000, 100000, by = 10000),
            labels = breaks / 1000000) +
  scale_y_continuous(name = 
            "Area within a five-minute biking distance\nof a public elementary (K-3) school\n(square km)",
            breaks = breaks <- seq(0, 800000, by = 100000),
            labels = breaks / 1000000) +
  theme_bw()